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A Multigrid Method for Two and 
Three-Dimensional Flows Behind Bluff 
Bodies in a Channel 


Abstract 


An investigation has been undertaken on the possibility of increasing speed of convergence 
of a existing computer code for the simulation of flow past bluff bodies in a channel. A 
multigrid method has been developed for the pressure-velocity iteration of the Marker 
and Cell ( MAC ) method used for the solution of full Navier-Stokes equations. For 
flows with complex configuration a remarkable reduction of computational time has been 
observed in a single V-cycle vis-a-vis single grid computations. The multigrid technique 
is based on the strategy that the low frequency errors of one grid can be solved on a 
coarser grid at the high frequency rate. The finer grid solution can then be corrected for 
the low frequency errors. The flow simulations using the multigrid scheme compare well 
with the single grid and available experimental results. 
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Chapter 1 


Introduction 


1.1 Background of the Problem 

A multigrid technique has been implemented to accelerate the convergence of the solu- 
tion of flow equations in a complex configuration. The flow in a channel with built-in 
square obstacle has been presented as a test problem for implementation of the multigrid 
technique. 

Fluid mechanics is an all pervasive subject. Application of fluid mechanics are 
numerous to list. In engineering, almost every equipment has some components that need 
to be cooled, heated, lubricated or supplied with energy. So, design of these components 
requires the knowledge of detailed velocity, temperature and pressure fields. Only detailed 
knowledge can lead to improved designs, low maintenance cost and reduced pollution 
etc.. Experiments can also give the detailed information about these parameters but the 
difficulties are with (i) cost, (ii) safety, (iii) feasibility and (iv) time. 

Predictive procedure is one of the solutions of eliminating these difficulties. But 
accuracy of a predictive procedure is very much dependent on the number of grids in 
the domain. If the number of grids are increased considerably, the computation becomes 
costly and obviously, the calculations take enormous time. In order to circumvent the 
reduced computational speed, the multigrid techniques are suggested. It is conjectured 
that the multigrid techniques will accelerate the computational speed well above the 
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computational time due to single-grid codes. 

In the present work, the flow in a channel with built-in obstacles is simulated by 
solving the full Navier-Stokes equations for incompressible viscous flows. The determi- 
nation of the dominating frequency in flow induced vibration with moderate Reynolds 
number is the main concern while designing big structures which are susceptible to os- 
cillations. The present study provides a fast flow-solver algorithm for determining the 
dominating frequency for such structures. 


1.2 Review of Literature 

Despite substantial progress in both hardware and software the computational fluid dy- 
namics calculations need higher computational speed. Highly accurate flow-physics sim- 
ulations clearly require enormous number of points in the computational domain which 
slows down the convergence. A question arises whether improvements in numerical tech- 
niques can be expected to show somewhat increasing convergence rate for a given accu- 
racy. 


The sequential scalar processors used in present-day hardware may be referred to 
as being single-instruction, single-data ( SISD ) architecture. The technological limits 
for performance enhancement of SISD are approaching rapidly. Vector processors use 
single-instruction, multiple data ( SIMD ) architecture to perform a single operation si- 
multaneously on a collection of operands. Parallel processors use multiple-instruction, 
multiple-data ( MIMD ) architecture to perform simultaneously several distinct opera- 
tions, with each operation accessing its own collection of operands. 

Together with the hardware improvement, one notable trend in fluid dynamics al- 
gorithm design has been the evolution of efficient algorithmic improvements, such as, 
domain decomposition and multigrid techniques. The multigrid method is indeed capa- 
ble of accelerating the speed of computation. The multigrid algorithms have been well 
documented by Brandt ( 1984 ) and Stuben and Trottenberg ( 1982 ). Multigridding is 
used to designate a process that carries a residual through a series of operations in order 
to accelerate the removal of error from the solution. The operations are related to the 
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grid size, which becomes increasingly coarse during one part of the iteration ( operation ) 
and increaisingly fine during the other. 

It is a common phenomenon in any type of relaxation scheme that after first five or 
six iterations the convergence rate slows down. The rate of convergence strongly depends 
on the mesh size. It becomes slower as the mesh is refined. At any stage of relaxation 
the error spectrum will contain a number of discrete frequencies, represented by Fourier 
components as shown in Fig. 1.1. 


k = 1 



Each frequency converges at a different rate. The high frequencies converge the 
fastest and the low frequencies the slowest. So, the convergence is held up by the longer 
wavelength or the low frequency components of the error spectrum. The nature of error 
elimination for different frequencies are shown in Figures 1.2(a) and 1.2(b). Each iter- 
ative procedure has different rates of convergence ( smoothing rates ) for each range of 
frequency. Out of the different methods of iterative ( relaxation ) procedures, Gauss- 
Seidel is the cheapest in workcount, but the convergence rate for low frequencies is slow. 


It is known from previous discussion that the low and high frequency demarcation 



Error 



23 50 75 100 

Iterations 


(a) 


Low freqneiacies: High frequencies: 



Before relaxation 



After relaxation 

reduction of the amplitude S ignificant reduction of the amplitude 

(b) 

Figure 1.2: Typical error smoothing behaviour of appropriate relaxation method 
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is a function of the grid. That is, the low frequencies on one fine grid can become high 
frequencies on a coarser grid. In multigrid technique, several grid levels are generated. 
Each grid level can be constructed by halving the fine grid or by doubling the coarse 
grid. It is easy to have grids of increasing fineness if one starts from coarsest grid. So, 
after generating all the grid levels, the errors in the solution are resolved efficiently by 
using these grids in a continuous cycling procedure. 

As a summary, it can be stated that the multigrid technique uses the following 
steps : 

1. Detecting the high and low frequencies. 

2. Resolving the high frequencies. 

3. Calculating the low frequency components. 

4. Transferring the residuals or defects to coarser grids. 

5. Solving the equations on coarser grids. 

6. Transferring the corrections back to finer grids. 

A large amount of literature ( Brandt, 1984; Stiiben and Trottenburg, 1982; Ghia ct 
al, 1982; Brockmeier et al, 1986; Vanka, 1986a ) exists on the theory and performance of 
the multigrids on model equations. However, despite its potential, a few investigations 
have been reported in which practically important viscous flows have been computed with 
the help of multigrid technique. A few researchers reported about the investigations of 
the driven cavity problem ( Ghia et al^ 1982; Agarwal, 1981; Philips and Schmidt, 1985 ) 
in which the stream function and vorticity equations ( Ghia et al , 1982 and Agarwal, 
1981 ) and the primitive variable equations ( Philips and Schmidt, 1985 ) are solved. 
Applications of multigrid technique to a sudden-expansion flow ( Philips et al , 1985 ) 
and a three dimensional duct flow ( Fuch and Zhao, 1984 ) have recently been reported. 
Multigrid solution of the Euler equations have been reported, by .Jameson and coworkers 
( 1983 ) and by Chima and .Johnson ( 1985 ). However, many practical viscous flows 
( such as those in combustors, turbomachinery, heat exchangers, etc.) are currently being 
solved by single-grid techniques. The application of the multigrid technique to accelerate 
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the convergence of these calculations remains to be explored. In recent works of Vanka 
( 1986a &: 1986c ), a calculation procedure based on a coupled simultaneous relaxation 
of the momentum and continuity equations was successfully used in association with 
the multigrid technique to analyse the complex flow fields in shear-driven square and 
cubic cavities. The relaxation procedure used in cavity problems is called Symmetrical 
Coupled Gauss-Seidel ( SCGS ) and by this procedure the equations are solved node by 
node in a coupled manner. At each node, the four ( or six ) velocities located on the 
cell faces and the pressure are simultaneously updated by the coupled solution of the 
relevant momentum and continuity equations. Such a coupled relaxation is advocated 
on the basis of earlier w'orks ( Vanka, 1985 k 19S6d; Galpin, 1985; Zedan and Schnei- 
der ,1985 ) that in internal flows, where the pressure field plays an important role in 
flow distribution, the rate of convergence is significantly improved when the momentum 
and continuity equations are relaxed simultaneously. In another investigation of Vanka 
( 1986b ), the performance of the above mentioned multigrid technique is assessed in four 
various complex configuration of recirculating flows important for industrial applications. 
The procedure is observed to converge rapidly to an acceptable accuracy in all the flow 
situations. 


1.3 Scope of the Work 


The objective of the present study has been to study the convergence properties of the 
multigrid strategy in complex flows, representative of those in many practical situations. 
The multigrid technique has been introduced only in the pressure-velocity ( continuity ) 
equation to get a faster convergence of correct pressure fields. This is the main condition 
to be fulfilled for obtaining divergence free flow fields. The multigrid algorithm has been 
tested through a test problem which has established single grid results. Single-grid results 
with respect to the vortex-shedding and wake-zone aerodynamics have been compared 
favourably with the multigrid results. 
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1.4 Organisation of the Thesis 


Chapter-1 of the thesis provides the introduction of the problem and a review of literature 
relevant to the basic multigrid theory with brief information about the scope of the 
work. The mathematical modelling of flow problem for numerical simulation has been 
presented in Chapter-2. Chapter-3 details the solution procedure to obtain the results. 
Implementation of the multigrid strategy has been discussed in Chapter-4. The numerical 
results has been presented in Chapter-5. Comparison of performance of the multigrid and 
single-grid has been accomplished in this chapter. The concluding remarks with scope of 
further research has been mentioned in Chapter-6. 



Chapter 2 


Mathematical Modelling of Flow 
Problem 


2.1 Introduction 


The oscillation that we generally encounter in chimney stacks and other vertical struc- 
tures in transverse flows is caused basically by vortex shedding. An initially smooth and 
steady flow across bluff bodies such as square cylinder, circular cylinder etc. may bring 
about the damaging deflection of the body in cases where the natural frequency of the 
obstacles are very close to the shedding frequency of the vortices. If the wake region 
behind the obstacle where Karman vortex street has been formed, is observed carefully 
it will be seen that the wake zone undulates like a flag from side to side. Due to these 
alternating deflections of the wake, periodicity is induced in the entire flow field behind 
the obstacle. As a consequence of this periodicity the forces on the bluff body become 
periodic in nature and Anally set the obstacle in vibration. If the excitation frequency 
of the force synchronizes with the natural frequency of the bluff body, the phenomenon 
of resonance will take place . Hence the study of unsteady flows across the bluff bodies 
have its relevance to design of structures, road vehicles, heat exchangers and wherever 
there is a possibility of flow induced vibration. 
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2.2 Effect of Blockage Ratio, Reynolds Number 
and Neighbouring walls on Vortex Shedding 
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Effect of various parameters on vortex shedding can be explained on the basis of the 
numerical investigation by Kiehm, I\Iitra and Fiebig ( 1986 ). They explained the solu- 
tions of the two and three dimensional Navier-Stokes and continuity equations for un- 
steady incompressible laminar flows behind a circular cylinder in a channel of rectangular 
cross-section at different Reynolds numbers and blockage ratios. Their investigations on 
two-dimensional flow show that the periodicity of the wake becomes apparent at higher 
Reynolds number with lower ratios of channel width H to cylinder diameter D, viz., at 
ReD = 53 for H/D = 5, and at Ren = 170 for H/D = 2. The channel walls show damp- 
ing effects on the periodicity and a parabolic flow profile evolves at the exit of a long 
channel. From their investigation, it can be concluded that neighbouring walls always 
provide a damping effect on the vortex shedding behind a bluff body. For a particular 
blockage ratio, the amplitude of periodicity increases with increase in Reynolds number. 
On the other hand, critical Reynolds number for vortex shedding increases with increase 
in blockage ratio. 


2.3 Governing Equations and Boundary 
Conditions 

In order to have quantitative information on the flow field in a channel with a built-in 
square obstacle and to get their dependence on various governing parameters, numerical 
simulation has been carried out. The geometry of the present problem is presented in 
Figures 2.1(a) and 2.1(b). 

The flow is simulated by solving unsteady Navier-Stokes equations for laminar, 
incompressible, two and three dimensional flows. The nondimensional equations for 
continuity and momentum for laminar incompressible flow may be expressed in their 
weak conservative forms as follows: 

at/ dV dW 

~ dX dV'^ dZ 


0 


( 2 . 1 ) 
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dU 

dip dUV dUW 

dP 

VPI 

(2.2) 

dr 

+ + -1- = 
dx dY dz 

'dX^ 

Re 

dv 

dVU dV^ dvw 

dP 

w 

(2.3) 

dr 

+ -h + = 

dX dY dZ 


Re 

dW 

dUW dVW dw^ 

dP 

VHV 

(2.4) 

dr ^ dX ^ dY ^ dZ 

dZ ^ 

Re 


In the above equations the velocities are nondimensionalized with the average in- 
coming velocity Uav at the inlet, all lengths with the channel height H and the pressure 
with pUav^- Reynolds number is denoted by Re. 


2.4 Boundary Conditions 


The boundary conditions of interest in the investigation are 


• Top and bottom plates: 


U = V = W = 0 : 


• side plates: 


U=V=W=0 ; 


• Channel Inlet: 


U = U{Y) , V = W = 0 ; 


• Channel Exit : 

There is no unique prescription for outflow conditions. The idea is to have such 
conditions which do not affect the flow in the upstream. The second or higher order 
derivatives of all the velocity components in the streamwise direction are put to 

d^U d^V dHV 
dX^ ~ dX^ ~ dX^ 


zero. 


(2.5) 
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Figure 2.2: Discretization of a tliree-dimensional domain 


• Obstacle : 

No-slip boundary conditions are used for the velocities on the obstacle. 


2.5 Grid System Used 

For single grid system, computational domain is divided into a set of rectangular cells 
( Fig-2.2 ) and staggered grid arrangement ( Fig.2.3 ) is used in which velocity components 
( vector quantities ) are defined at the centre of the cell faces to which they are normal. 
The pressure and the divergence ( scalar quantities ) are defined at the cell centre. The 
advantages of the staggered grids are : 


• pressure difference between two adjacent cells is the driving force for the velocity 
components located between the interface of these cells. 

• transport rates across the faces of the control volume can be computed without 
interpolation of velocity components. 



Figure 2.3; Three-dimensional staggered grid showing the locations of the discretized 
variables 

2.6 Numerical Boundary Conditions 

We have discussed the boundary conditions of symbolic form. It is necessary to dis- 
cuss the implementation technique of the boundary condition in the solution code. The 
boundary conditions are imposed by setting the appropriate velocities in the fictitious 
( or imaginary ) cells surrounding the physical domain. 

2.6.1 Boundary Conditions for Confining Walls 

The governing differential equations are elliptic in space and parabolic in time. We need 
the boundary conditions for all the confining surfaces. Since the confining walls at the 
top and bottom are rigid for the present computation, both the normal and the tangential 
components of velocity on the wall must be zero. With reference to Fig. 2. 4, we can write 
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( for the bottom boundary of the computational domain ) 


= —Ui^2,k 

Km, A- = 0 

VKm,/,- = -1'K,2.x- 


for i 
and k 


2 to ire 
2 to kre 


( 2 . 6 ) 


Similarly, for the no-slip right wall, the boundary conditions become 


K,i.i 


-Uid,2 

-Kj.2 

0 


for i = 2 to ire 
and j = 2 to jre 


(2.7) 


The front plane ( inlet plane ) is to be provided with the inflow boundary condition. 
The normal velocity components are set to zero and a uniform or parabolic axial velocity 
may be deployed. Hence, with reference to Fig. 2. 4, we can use 


If 1 . j,fc 

K,m- 

,j,k 


-K2,M- 

1.0 


1.5 


1 - 


jm j 

Jm 


2' 



j 

k 


2 to jre 
2 to kre 


( 2 . 8 ) 


where jm is the horizontal mid-plane. 


Continuative or outflow boundaries always pose a problem for low-speed calcula- 
tions, because whatever prescription is chosen it can effect the entire flow field upstream. 
We should follow a prescription that allows fluid to flow out of the computational do- 
main with a minimum of upstream influence. In our present computation, the second 
derivative of the dependent variables with respect to flow direction are set to zero to 
ensure smooth transition through the outlet boundary. These can be implemented in the 
following way: 


Ui+l,j,k 


- 1 

for 

j = 2 to 

jre 


K+1 J,X: 

= 2K,,m 

- K-n,.. 

} and 

k = 2 to 

kre 

(2.9) 


= 2W.-JM- 

- J 

at 

i = iim — 

2 ' 



The governing equations and the boundary conditions discussed so far are described for 
a three dimensional problem. The two dimensional problem is a subset of the three 
dimensional problem and the implementation is straightforward. 
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J=jim 
j = jre 


Pii 


i = iim 


jiiiPSPSPSPSPSP^iiil 


k=kim^“^‘'® 8z i«/k=l j 


Right side wall 


Out flow boundary 


Bottom boundary 


Front inflow plane 


Figure 2.4. Boundary conditions and fictitious boundary cells 

2.6.2 Boundary Conditions for the Obstacle 

In order to implement the kinematic boundary conditions on the obstacle, the velocity 
points which are in the domain of the obstacle are set equal to zero. The implementation 
IS explained with the help of a three dimensional flow simulation. Application for a 
two dimensional case is straight forward. Implementation of the no-slip condition for V- 
and W-components of velocity along the transverse direction, U-component of velocity 

along the direction of flow need some manipulations. The boundary conditions are as 
the following ( see Fig. 2. 5 ) : 

for i = ia — I to ih at k = ka 
and = ~L/{j^k+i for i = ia — I to ib at k = kb 

similarly, for V-component of velocity, we have, 

= —K’j.A;-! for i = ia to ib at k = ka 


( 2 . 10 ) 

( 2 . 11 ) 


( 2 . 12 ) 
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Kj.fc — fo7' 

i = ia to 

ib 

at 

k = kb 

(2.13) 

Vi,j,k = for 

Ic — - IcCL to 

kb 

at 

i = ia 

(2.14) 

and Vij^k = - Vi+i,j,k for 

k = ka to 

kb 

at 

i = ib 

(2.15) 

finally, for W-component of velocity, one can have 





I'Kj.A: = for 

k = ka — \ 

to 

kb 

at i = ia 

(2.16) 

(ind ^'V{jk J or 

k = ka — 1 

to 

kb 

at i = ib 

(2.17) 

For all the above cases, j varies from 2 to 

jre 







Figure 2.5: Discretization of the three-dimensional square obstacle 


2.7 Closure 


In the present numerical investigations, the input parameters are Reynolds number, 
nondimensional length of the channel, inlet velocity profile and the nondimensional ob- 
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stacle height. For two-dimensional problem, the nondimensional length of the channel 
is 4.875. The obstacle is placed at a nondimensional distance of 1.375 from inlet. In 
three-dimensional problem, the nondimensional breadth and height of the channel is 1. 
The aspect ratio of the obstacle is 1. The nondimensional length of the channel is 3.333. 
The obstacle is placed at a location whose nondimensional distance from inlet is also 1. 



Chapter 3 


Solution Procedure 


3.1 Solution Algorithm 

The nondimensional continuity and momentum equations written in terms of primitive 
variables in weak conservative forms ( equations(2.1) through (2.4) ), are solved by using a 
modified version of Marker And ^ell (MAC) algorithm. The original version of the MAC 
method due to Harlow and Welch ( 1965 ) has been modified by Hirt and Cook ( 1972 ). 
In the original MAC method, the pressure field is determined by directly solving Poisson 
equation for pressure, whereas, in the modified version of MAC algorithm pressure values 
are calculated by a pressure-velocity iteration procedure. A related technique developed 
by Chorin ( 196T) involved a simultaneous iteration on pressure and velocity components. 
Viecelli ( 1971 ) has shown that the tw'o methods as applied to the MAC algorithm are 
equivalent. 


3.2 Details of the Computational scheme 

The MAC is a semi-implicit scheme for solving full Navier-Stokes equations where the 
advancements of velocity components in the time direction are obtained explicitly by 
calculating accelerations due to convection, diffusion and pressure gradient. Thus, after 
obtaining provisional velocity field for a time step, the continuity equation is solved 
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implicitly to correct the provisional velocity field. Hence MAC is a semi-implicit scheme. 

From the guessed velocity and pressure fields, the corrected velocity and pressure 
fields are obtained by pressure- velocity iteration through the continuity equation. Con- 
vergence of this iteration process ensures a divergence-free velocity field for the initial 
time step. Now, the corrected pressure-velocity fields are used to calculate the velocity 
field for the next time step by solving the Navier-Stokes equations. The time step 5t for 
the adv'ancement of the velocities is governed by the stability criteria which are discussed 
in a subsequent section. However, this explicit time advancement may not lead to a 
velocity field with zero mass divergence in each cell. The satisfaction of the continuity 
equation is achieved by adjusting the pressures as well as velocities in each cell through 
an iterative process, which as mentioned earlier, is equivalent to solution of a Poisson 
equation for pressure. Details of pressure- velocity iteration procedure is discussed in the 
following section. The corrected pressure and velocity fields are now used to evaluate 
the provisional velocities for the next time step and so on. Thus, starting from initial 
guess fields, the solution advances in time direction until a steady or periodic solution is 
obtained. 


3.3 Discretization of Governing Differential 
Equation 


Because of staggered grid arrangement ( Fig.2.3 ) the velocities are not defined at the 
nodal points. Whenever required, the velocities at the nodal points are calculated 
by interpolation of neighbouring velocities. As an example, we can write = 

0.5((/,_ij,jb -f but where a product or square of such quality appears, it is cus- 

tomary to interpolate first and then the product to be taken. 


For any cell, the discretization of the continuity equation is given by 




8X 


T ;:77 1 


SY 


SZ 


(3.1) 


The discretization of the temporal derivatives in equations (2.1) to (2.4) is done by 
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a forward differencing 


au 

dr St 


(3.2) 


The convective terms of the momentum equations are discretized by a weighted av- 
erage scheme which combines the upwind and central diiferencing to achieve the stability 
of upwind and formal accuracy of the central differencing ( Hirt, Nichols and Romero, 
1975 ). The discretization of one of the convective terms is shown below; 


d{UV) 

dY 


+ + K + l,j,A:)l — Ui J + l,k) 


(3.3) 


The factor Ofp determines the relative weightage of the central and upwind differenc- 
ing. ttp = 0 gives simple a central differencing and Op = 1 leads to a second upwind 
differencing. The value of this factor is between 0.2 and 0.3. 


3.4 Pressure- Velocity Iteration 


In the original MAC method the Poisson equation was solved for pressure. The mod- 
ification due to Hirt and Cook ( 1972 ) enables one to solve the Poisson equation for 
pressure in such a way that the velocity and pressure fields are corrected in an iterative 
procedure. 

The velocities are advanced explicitly for the next time step from the momentum 
equations if the velocities and pressures are known at any time step. But these velocities 
do not necessarily describe a meaningful solution if there is a non-zero divergence in any 
cell. The pressures distribution are correct if it is accompanied by a zero divergence 
in each cell. The presures have to be corrected so that there is no accumulation or 
annihilation of mass anywhere in the domain. Corrections to pressures are applied in an 
iterative manner as described by Hirt and Cook ( 1972 ). 

The method of pressure-velocity iteration is as follows. The relationship between 
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the explicitly advanced velocity components at (n + l)th time level and those at the 
previous nth time level can be written as 




= + 



+ 5r[RES!DU]l^^^ 


(3.4) 


where [RESIDU]lj^ j, is the discretized form of 

dU^ dUV dUW V^U' 

~~dx~'w ^^~Rr 

at the {i^j,k) cell evaluated with the velocity values at the nth time step. 


On the other hand, the corrected velocity components, which are still unknown, are 
-related to the correct presure ( unknown ) in the manner 


+ (^) + sriREsiDUYi,, 


from equations (3.4) and (3.6) one can write 


(3.6) 




(3.7) 


where Pr, = , 


Neither nor are explicitly known at this stage so that one can be calcu- 
lated with the help of the other. Calculations are done in an iterative cycle and we can 
write 



In the similar manner the other velocities can be written as 


Corrected 




Estimated 

- 


KS + 


Correction 

(If) 

1^) r - /='..i+u) 


(3.8) 


(3.9) 

(3.10) 

(3.11) 
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+ [ 7^ \{P'i.,*-P'. ■,.<:») 


— '•^’S-l - ' (P'i.i/:-P'i.,J^-^) 



( 3 , 12 ) 


( 3 . 13 ) 


( 3 . 14 ) 


The correction is done through continuity equation. Plugging in the relationship 
into the continuity equation we get 

SX SY SZ 

SX SY SZ 


— 8r 


— Sr 


SX^ 

■P/j+l.A- ~ ‘^Pjj.k + -P/j-l.fc 

6Y^ 


At this stage, the pressure corrections in the neighbouring cells are neglected. There- 
fore, we get the velocity corrections in the form. 









yn+i ^ all 





iP'iJJc) 


{P'i.iJ:) 








( 3 . 16 ) 


( 3 . 17 ) 


( 3 . 18 ) 

( 3 . 19 ) 


( 3 . 20 ) 


( 3 . 21 ) 
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The continuity equation takes the following form after neglecting the pressures cor- 
rections in the neighbouring cells, 

- <77-, , irS - vrgiL, 

<5A' ^ 5V" ^ 6Z 

5X 5Y 5Z 


or, 


or, 


8t 


+- ( ^ ) M ^ + 2 ( ^ ) {pu:] 


St 


SY^ 


5Z^ 


0 Du,k + Pij,k 2(5t I 


Pi,j,k - -^'oi Di.j,k)/ 


\r f 1 1 J 

TVS + + 


[SX'^ SY'^ SZ'^il 


(3.22) 

(3.23) 

(3.24) 


where ojo is an over relaxation factor which is introduced to accelerate the pressure 
correction process. Usually a value of 1.7 is used. After calculating velocities in 

each cell are corrected according to the equations (3.16) to (3.21) and the pressure in 
each cell is adjusted as 




Pi'j.k + PUk 


(3.25) 


This process is continued until the velocity divergence in each cell vanishes. If the 
velocity boundary conditions are correct and a divergence-free converged velocity field is 
obtained, eventually correct pressures will be evolved in all the cells including the cells 
on the boundary. This feature of modified MAC method has been discussed in more 
details by Peyret and Taylor ( 1983 ). However, it was also shown by Brandt, Dendy and 
Ruppel that the aforesaid pressure-velocity iteration procedure is equivalent to solution 
of Poisson Equation for pressure. 


3.5 Stability Criteria 


For a given mesh the choice of the time step is determined through stability analysis which 
has to take care of two conditions. First, fluid should not allowed to cross more than 
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one cell in one time step . This restriction is derived from the Courant-Friedrichs-Lewy 
(CFL) condition given by 

, . {8X SY 8Z\ , ^ 

i,, (3.26) 

where the minimum is with respect to every cell in the domain. Typically 8t is chosen 
equal to one-third to two-third of the minimum cell transient time. 


Secondly, when a nonzero value of kinematic viscosity is used, momentum must not 
diffuse more than approximately one cell in one time step. A linear stability analysis 
shows that the restrictions on grid Fourier number will yield 


1 {8X‘^){8Y'^){8Z'^) 


(3.27) 


Finally, the minimum of the above two time-steps is chosen for the computation. 

The term ctp in the discretization of the convective terms of the Navier-Stoices equa- 
tions gives the desired amount of upstream ( donor cell ) differencing. When oc^ = 0, 
a space-centred differencing is obtained and Op = 1 leads to full upstream or donor cell 
form. In general, should be chosen slightly larger than the maximum value of 


U8t 


8X 


or 


V8r 


8Y 


or 


W8t 


8Z 


(3.28) 


occuring in the entire domain. In other words, 


1 > cXp > max ■ 


U8t 


8X 


V8t 


8Y 


W8t 


8Z 


(3.29) 




Chapter 4 


Multigrid Implementation 


4.1 Objective 


It has been experienced that most of the computational work of the MAC algorithm 
is due to the pressure-velocity iteration. The purpose of our investigation is to reduce 
computational work by the application of a multigrid techniques to the pressure-velocity 
iteration. In order to explain the computational strategy of the multigrid technique, we 
shall start with the basic discretized equations which have been explained in Chapter-3. 
However, we shall use the symbolic forms which are more convenient for explaining of 
multigrid strategy. The difference equations approximating the momentum equations are: 



- 

+ 

(SrlSX)(P2S.i.k- PS) = 



VIS 

- 


(i^liy)(SS.k-PS) = 


> (4.1) 

vvm 

- 

+ 

(Srl6Z)(PSS-PS) = 

J 



where the superscripts n stands for the time step ; a, b, c represent the combined 
convective-diffusive terms. The continuity equation is discretised as follows: 


= jx + SY + Jz “ 


The computational sequences of equations (4.1) and (4.2) is shown in Fig.4.1. 
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Initial pressure and 
velocity at time X 


Provisional velocity field 
from the N-S equation 


Pressure - velocity iteration 
through continuity equation 
by a number of SOR 
sweeps 


r= At 


Time to stop 


Print 


Stop 


Figure 4.1. Computational sequence of the A'lAC algorithm 
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The treatment of the continuity equation together with the pressure-velocity itera- 
tion in MAC method is equivalent to a solution of Poisson-equation for pressure which 
is given by 


P'{X) 


-u.<oD{X) 


2St 


1 1 
SX-2 ^ 


+ 


5Z'^ 


( 4 . 3 ) 


The Gauss-Seidel sweep is performed to reduce D in each cell of the computational 

grid. 


During the initial stage of computations, a large number of relaxation sweeps for 
each time step is necessary in order to reduce the divergence ( error ) for some initial 
value Dj below a given upper bound Dp{^ 10““*). 

These errors reduced by such relaxation methods may roughly be divided into low 
and high frequency components corresponding to meshsize h. The convergence usually 
slows down after the high frequency components being smoothened out. Convergence 
properties for the reduction of a low frequency error component corresponding to h can 
be improved by treating this component on a coarser grid with meshsize 2h where it 
appears as high frequency component corresponding to 2h. 


4.2 Grid Arrangement 

The multigrid technique uses several grid levels in cyclic order to accelerate the conver- 
gence of the solution. But the generation of grids needs special care such that it does not 
affect the boundary conditions. Generally, the calculation is started with the coarsest 
grid. Then this coarsest grid is splitted into several levels by doubling the number of 
grids ( or halving the mesh size ). The finest grid is used for solving the basic difference 
equations. Calculation can also be started from the finest grid depending on the type of 
cycling scheme and equations to be solved. A typical two-dimensional grid arrangement 
is shown in Fig. 4. 2 for three different grid levels. An extension to three dimensional 
geometry is straight-forward. 
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4.3 Multigrid Strategy 


In multigrid method error components of different frequencies are reduced on a number of 
grids having successive coarseness, Brandt et al ( 1980 ) suggested the multigrid technique 
of equation (4.2) can be carried out by cori'ection scheme (CS) mode because of the linear 
character of the differential operator of the continuity equation. 

Several cycling schemes for fine to coarse grids are possible. In one strategy V-cycles 
are used ( Fig.4.3 ) to compute equation (4.3) on a fine grid and a number of correction 
equations (4.5) and (4.6) on successive coarser grids. 


A V-cycle starts from the fine grid and performs a previously specified number of 
iterations on the fine grid and suljsequent coarse grids. RESTRICTION is the process of 
obtaining a coarse grid solution from a finer grid. A V-cycle ends at the finest grid by per- 
forming iteration after corrections are prolongated on each grid level. PROLONGATION 
refers to extrapolation of coarse grid values to finer grids. The symbols which we shall 
follow in the subsequent sections are as follows: 


P 

Ui 

h{6X,6Y,5Z) 
and X 


N ondimensional pressure field 
{U, V, W) = Velocity field 
Mesh size 
Cell center 


4.3.1 Restriction 

To start with the Gauss-Seidel relaxation sweeps are done with the correction equation 
(4.3) on finest grid. At this stage we can write 

Di = div Ui on finest grid (/ = m) (4-4) 

where, div is a divergence operator 

Subsequently the Gauss-Seidel relaxation sweeps are accomplished on the coarser 
grids (/ < m). The operation can be stated as 
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Figure 4.3: A V - cycle ( CS mode ) 



Figure 4.4: R.e.striction of residuals ( defects ) to coarse grid 
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Di = div Ui + di starting with Ui = 0 ] Pi = 0 (4-5) 

where, di = sum of restricted defects ( divergence ) of grids from I to (/—I), / being < m. 

After a specified number of iterations at a particular grid level, next coarser grid 
will be taken up and the divergence ( defect ) is restricted for the next grid. The general 
algorithm is: 

m 

Computation of defect : di = '^Di (4.6) 

i 

Restriction of defect[divergence) : di = ll~^di (4-7) 

where, is restriction operator which is explained in Fig.4.4. Numerically the operator 
can be expressed as 

= (4.8) 

^ n 

The process of taking up next level of coarser grid continues till the coarsest grid is 
reached. RELAXATION until convergence is carried out only on the coarsest grid (1=1), 
where the calculation is cheap because of the number of grid points. Experience shows 
that a number of sweeps <3) are sufficient at each grid level. However, such specified 
number of sweeps at each grid level may require multiple V-cycIes. A typical value of 
number of V-cycles has been reported as 8 ( see Brockmeir et al, 1986 ). In order to 
reduce the number of cycles, an adaptive cycling scheme is suggested ( Vanka, 1986b). In 
adaptive cycling scheme, a fixed error reduction is desired rather than a fixed number of 
iterations. 

Let be the error on grid I at kth. iteration, then switching over to (/ — 1) grid is 
recommended if is achieved. We have used the adaptive cycling scheme 

and the value of has been chosen as 0.6. Iterations are performed on grid I until the 
above criterion is met. 

For two dimensional problem the restriction can be explained with the help of 
Fig. 4. 5. Let {ic,jc) and {if,jf) are the coarse and fine grid indices respectively. One 
can write 
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if = 2 X ic — 1 
and jf = 2 X jc — 1 

The restriction operator explained by equation(4.8) can be written as 

d^'iicjc) = 0.25 X [d^{if,jf) +d^{ifjf - 1) + 

d^if - l,i/) + d^{if - IJf - 1)] (4.9) 

However, while switching over to the next coarser grid, the velocities and pressure at 
each grid level are stored in a separate array. 

4.3.2 Prolongation 

As it has been stated earlier, after having achieved convergence in coarsest grid, the 
values of pressures are prolongated to the next finer grid. Weighted interpolation is 
used to obtain finer grid pressures corrections out of coarse grid pressures ( Hackbuch 
and Trottenberg, 1982 ). For simplicity this interpolation is shown schematically for a 
two dimensional geometry in Fig.4.6. Most commonly used interpolation is trilinear. A 
trilinear technique considers a linear variation in each spatial direction. For a uniform 
two dimensional grid, one can write 

/a = (l/16)(9i}+3C + 3A + 5) ' 

/2 = {1/16){9C + 3B + 3D + A) 

f (4-10) 

/3 = (1/16)(9A + 3S + 3P + C) 

U = (1/16)(95 + 3A + 3C + P) , 

In terms of operators, as it was done in case of restriction, the equations (4.10) can 
be expressed as 

l^P'm = + ^P'b + ^P'b + P'c)i-i 


( 4 . 11 ) 
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The corrections of the finer grid pressure and velocity approximations are then 
carried out in the following way 

AP/(X) = 

Pl(X) ^ Pl(X) + APl(X) (4.12) 

U,(X + h/2) ^ Ui(X + h/2) + (AT/h)[^\Pl(X) - APl(X + h)] ^ 

After prolongation, only one sweep ( 1^2 = 1) is good enough for each grid level. 



■ 27 9 9 9 3 3 3 1 ' 

9 27 3 3 1 9 9 3 

9 3 27 3 9 1 9 3 

ri _ 9 3 3 27 9 9 1 3 

~ 64 3 1 9 9 27 3 3 9 

3 9 1 9 3 27 3 9 

3 9 9 1 3 3 27 9 

_ 1 3 3 3 9 9 9 27 _ 

As it can be seen from Fig. 4. 6, pressure corrections in the cells on the boundary can 
not be calculated out of coarse grid pressure corrections ( P/_j ). The pressure corrections 
at the boundary cells are obtained from the inner pressure correction terms (A P/s) by 
linear extrapolation normal to the boundary. 
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4.3.3 Velocity Boundary Conditions for Multigrid 


It has already been discussed that the velocity boundary conditions are needed to solve 
equation(4.3). The boundary conditions for pressure are not needed. 




G for I — m 
0 for I < m 


(4.14) 


and G = 5„(.T,T/,z) 

3v{x,y,z) i (4.15) 

9 w{x,y,z) ^ 

where and g^u are the boundary conditions for the integration domain and the 

subscript B stands for the boundary. 


4.4 Closure 

The important aspects of a multigrid computer program are the RESTRICTION, 
RELAXATION and PROLONGATION. These must be carefully coded and checked out. 
The total storage necessary ( excluding boundaries ) is the sum for all grids. For a two 
dimensional domain, this is 

771 n 771 n 
mxn + -x- + -x- + ... 

11 4m 

= m X n (1 + - + — + •••) = — X n 
4 16 3 

where m x n is the storage on the finest grid. Extension for a three dimensional domain 
is straight forward. 



Chapter 5 


Results and Discussions 


5.1 Flow Simulation 

The numerical investigations of the structure of confined wakes behind a square cylin- 
der for two-dimensional and three-dimensional cases are presented. The numerical flow 
visualization provides the details of the vortex-shedding phenomenon. 

Experimental and numerical studies have shown that the vortex shedding behind 
the cylinder introduces periodicity at moderate Reynolds numbers in the flow field. The 
nondimensional Strouhal numl^er which varies with Reynolds number as well as with the 
blockage ratio characterizes the unsteady periodicity of the wake. The damping of the 
periodicity in the downstream of the flow may be attributed to the influence of sidewalls 
on the flow structure. The present study has analyzed the fundamental flow structure 
and the nondimensional relationships characterizing the vortex-shedding frequency. 

5.1.1 Two-Dimensional Flow in a Channel With a Built-in 
Square Cylinder 

For the computation of the present problem, MAC algorithm with 314 x 66 grid has been 
used. Uniform grids are deployed throughout the calculation domain. Computations 
have been carried out in a channel of length Lj H — 4.875. The geometrical centre of the 
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— =Uav Res = 40 B/H = 0.25 



Figure 5.1: Attached vortices behind the obstacle at Res = 40 

square cylinder is located at a distance Xr = 1.375 from the inlet. The channel and the 
cylinder axes are aligned ( Fig. 2. 1(a) ). The aspect ratio of the square cylinder is one. 
The simulation is carried out with two different blockage ratios {BjH), namely, 0.125 
and 0.25. The Reynolds number based on the cylinder height, Rcr is used as the input 
parameter. 

The structure of the wake and its functional relationship with Reynolds number are 
studied by varying the Reynolds number for the same blockage and aspect ratios. In 
Fig.5.1, a steady solution for BjH = 0.25 and Reynolds number Rcr = 40 is shown. 
The recirculating wake is extended to nearly twice the obstacle width in the downstream 
of the obstacle. Figure 5.2 also shows a symmetric wake for a Reynolds number of 60. 
Comparison of the Figures 5.1 and 5.2 shows that the extent of the recirculating wake 
increases with the increase in Reynolds number. The flow field is steady in both the 
cases. The flow with Reynolds number of 60 does not indicate any vortex shedding. But 
at Rcr = 80, the wake loses its original symmetry and the flow becomes periodic in the 
near wake ( Fig. 5. 3(a) ). However, the periodicity is suppressed in the downstream of the 
square cylinder by the channel walls and the velocity profile at the exit of the channel 
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— Rcb = 60 B/H = 0.25 



Figure 5.2: Attached vortices behind the obstacle at Res = 60 

tends towards the steady parabolic ( not shown in Fig. 5. 3(a) ). For better understanding, 
a magnified view of the vortices behind the obstacle is shown in the Fig.5.3(b). Okajnna 
( 1982 ) observes the induction of the periodicity in the wake behind a rectangular cylin- 
der in an infinite medium ( blockage ratio being negligibly small ) at Rts = 70. Davis, 
Moore and Purtell ( 1984 ) finds periodicity in computations with the same geometrical 
configuration as the present one, at Res « 100. In another numerical investigation, 
Mukhopadyay d al ( 1992 ) observes the critical Reynolds number describing onset of 
periodicity as 85. Experiments due to Shair et al ( 1963 ) reveals the critical Reynolds 
number as 130 for the wake to become periodic in a channel with a circular cylinder 
in which the blockage ratio is 0.33 and the Reynolds number is defined in terms of a 
maximum velocity which would exist at the same location as that of the the centre of the 
cylinder under flow conditions identical with those of the experiment but m absence of 
the cylinder. It becomes clear that their critical Reynolds number describing the onset 
of periodicity, based on a definition as ours, would be 87. This is because of the fact that 
the average velocity is two-thirds of the maximum velocity. The critical Reynolds number 
obtained by present computation is well within the range of the Reynolds numbers ob- 
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- =Uav 


Res = 250 


B/H = 0.25 





Res = 250 


B/H = 0.25 



(b) 


Fieure 5.4: Velocity vectors in the channel at 
( Lb = 250 ) (a) r = 5.0823 (b) r = 7.4676 


different nondimensional time, 
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Figure 5.5: Streamline cros,sing the cylinder in the channel, ( Res = 375 ) 
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Obs. 

No. 

Res — 

pR av R / ^ 

B/H 

Strouhal Number, S 

From time- 
evolution of 
lift 

coefficient 

From signal 
traces of 
fluctuating 
velocity 

1 

80 

0.125 

0.1836 

0.1806 

2 

100 

0.125 

0.1876 

0.1829 

3 

200 

0.125 

0.1840 

0.1865 

4 

300 

0.125 

0.1819 

0.1838 

5 

400 

0.125 

0.1817 

0.1827 

6 

500 

0.125 

0.1800 

0.1808 

7 

600 

0.125 

0.1784 

0.1790 

8 

700 

0.125 

0.1773 

0.1779 

9 

800 

0.125 

0.1765 

0.1768 

10 

900 

0.125 

0.1749 

0.1752 


Table 5.1: Variation of Strouhal number with Reynolds number 


tained by other researchers. Finally, it is worthmentioning that the velocity vector-plots 
( Figures 5.1, 5.2 and 5.3 ) are not merely qualitative since the representative scales of 
the average channel velocity Uav are shown in the velocity vector-plots and these can be 
used for quantitative analysis of the results. 

With the increase in Reynolds number beyond the critical one, the von Karman vor- 
tex street is formed and alternate shedding of vortices into the stream becomes prominent 
( Figures 5.4(a) and 5.4(b) ). Figure 5.5 reveals that for a Reynolds number of 375, the 
flow separates at the leading edge and does not reattach. Needless to say that the entire 
wake zone undulates like a flag for this Reynolds number. At a low Reynolds number, 
there is a steady reattachment behind the leading edge and the flow finally separates at 
the trailing edge which results in symmetrical vortices. But at somewhat highef Reynolds 
number, after the separation at the leading edge, the flow reattaches either the upper or 
lower surface of the obstacle during alternate shedding of the vortices into the stream. 
However, in the range of Reynolds number of 375 the flow remains detached behind the 
leading edge. 
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Frequency (/) of vortex shedding lias been determined by a technique which borrows 
some idea from the measurements of Okajima ( 1982 ). The normal component of velocity 
( V-velocity ) was recorded at a location which is 3B behind the cylinder at the central 
axis of the channel. In case of steady flow, the flow smoothly divides and reunites 
around the cylinder, the consequence of which the normal component at that particular 
location will be zero. But at higher Reynolds number ( above the above mentioned 
critical Reynolds number ), the normal velocity keeps fluctuating at the same location. 
The frequency can be determined from the oscillating signal of the velocity component 
stored over a fairly long duration of time. Figure 5.6 shows the time evolution of the 
lift coefficient for two different Reynolds numbers. The vortex shedding frequency can 
also be determined from the time evolution plot of the lift coefficient distribution. The 
time period T can be calculated computationally by monitoring the non-dimensional 
time when the lift coefficient is just crossing the meanvalue. The difference between two 
such alternative time values ( as shown in Fig. 5. 6(a) and (b) ) gives the time period T. 
Once the time period T is known, the corresponding frequency /{= 1 /T) and the Strouhal 
number 5(= fBfUav) can also be calculated. It seems from the Figures 5.6(a) and 5.6(b) 
that the amplitude of the fluctuating coefficients increases with the increase in Reynolds 
number. Table 5.1 shows the Strouhal number calculations from two methods, namely 
the signal traces of the fluctuating components and the time-evolution of lift coefficients. 
The presented data shows agreement between the two values. Table 5.2 compares the 
Strouhal irumber variation with the experimental results of Okajima ( 1982 ). The present 
calculation gives somewhat higher Strouhal numbers for the entire range of Reynolds 
numbers. The discrepancy between the experimental results ( Okajima, 1982 ) and the 
present results is due to the fact that the experiments were done for a negligibly small 
blockage ratio. The finite blockage ratio {BfH = 0.125) might have brought about some 
changes in Strouhal number because of the influence of side walls. 

5.1.2 Three-Dimensional Flow in a Channel With Built-in 
Square Cylinder 


In three-dimensional flow simulation, 82 x 26 x 26 grid has been used at the finest level. 
Computations have been done in a channel, of length Lj H = 3.333 with a built-in square 



Cl( Lift Coefficient ) 
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Obs. 

No. 

Res = 

pUav^ / 

B/H 

Strouhal Number, S 

Experimental 
results in 
infinite 
medium 

( Okajima, 1982 ) 

( Upper limit ) 

Numerical 
prediction in 
present 
computation 
(From signal 
traces of 
fluctuating 
velocity ) 

1 ^ 

SO 

0.125 

0.130 

0.181 

2 

100 

0.125 

0.148 

0.183 

3 

200 

0.125 

0.150 

0.187 

4 

300 

0.125 

0.148 

0.184 

5 

400 

0.125 

0.142 

0.183 

6 

500 

0.125 

0.139 

0.181 

7 

600 

0.125 

0.135 

0.179 

8 

700 

0.125 

0.133 

0.178 

9 

800 

0.125 

0.130 

0.177 

10 

900 

0.125 

0.128 

0.175 


Table 5.2: Comparison of present computation with experimental ( Okajima, 1982 ) 
results 


cylinder of cross sectional area A y. B. In this problem, Ajli = BjH = 0.33 has been 
used. The obstacle is placed with its geometrical centre at a nondimensional distance 
Xr = 1 from the inlet. The channel and obstacle axes are aligned. The height of the 
obstacle is same as that of the channel. Only one blockage ratio, iiiz, 0.33 is used for this 
simulation. 

In Fig.5.7(a), the particle traces in a channel for a low Reynolds number flow Res = 
150 are shown. The plotted lines are the pathlines showing the way of fluid particles 
through the channel beginning at various starting points. This plot can be depicted 
as the numerical flow visualization. Under steady state, these lines are streamlines. 
However, in Fig.5.7(a), the structure of wake in the immediate neighbourhood of the 
horizontal symmetric plane is nearly the same as two-dimensional flows. 

However, in another view ( Fig.5.7(b) ) between the plane of symmetry and the 
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bottom wall, a screw-like helical structure of flow is observed. As such, Fig.5.7(b) explains 
the formation of a horseshoe vortex system. As fluid approaches the stagnation line of a 
square obstacle in cross flow, it slows down and its pressure increases. The smaller velocity 
in the boundary layer on the flat bottom plate which supports the square cylinder leads 
to a smaller pressure increase. Thus the induced pressure gradient causes a flow towards 
the bottom wall which interacts with the main-stream. The fluid rolls up forming vortices 
which are finally swept around the obstacle base, and are carried downstream. 


5.2 Computational Aspects of the Multigrid 
Technique 


In this section, the performance of the multigrid algorithm for both the two-dimensional 
and the three-dimensional flow situation with built-in square obstacle is presented. In 
order to assess the behaviour of the algorithm critically, several calculations with variation 
in the Reynolds number have been made for each of the flow situation. Since the multigrid 
strategy has been applied only to pressure-velocity equation, the results of the single- 
grid and multigrid should be the same. It has been observed from the velocity-vector 
plots that both the results are in good agreement with respect to their quantitative and 
qualitative measures. It is well known that the relaxation factor has its significant role 
in the convergence of solution. In single-grid calculation this factor is chosen as 1.8. But 
for the multigrid computation, special care has to be taken in choosing the relaxation 
factor. Successive over-relaxation has not been used as the smoothing operator. The 
over-relaxation destroys the high-frequency smoothing that is so crucial for the multigrid, 
method. In our present calculation a fixed relaxation factor is used equal to one. 


5.2.1 Two-Dimensional Case 

For the computation with the multigrid technique, grid size of 314 x 66 has been used 
as the finest grid. The number of grid level is resti'icted by obstacle size. However, we 
have computed for three number of grid levels with the boundary conditions in every 
grid levels. In Chapter-4, we have already mentioned about the intei'polation of values 
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Rcb 

~ pUayR ! p 


CPU Times 

Hr : Min : Sec 

Workxjm- amits 
to get CO ■merged 

steady s^;;;;;^lution 

Amount of 
work in 
relation to 
single grid 

25 

Multigrid 

1:10:46.4 

lOO 3_ 

11.98 % 

Single-grid 

2:05:17.9 ' 

836 

40 

Multigrid 

1:18:31.0 

111 

12.45 % 

Single-grid 

2:13:22.3 


50 

Multigrid 

2:10:31.2 


17.76 % 

Single- grid 

3:17:45.0 

115:^06 


Table 5.3: Comparison of CPU-times and Workunits bet\v^.^xi the multigrid and single- 
grid for two-dimensional problem 

near the boundary where the general prolongation rule doe:^ iiot hold good. But in our 
present calculation, we have made the prolongation routines such a way that it does 
not require any linear extrapolation of inner values to get boundary values. 

The present computation has been performed on a 0000/735 series computer. 

The system hcis 144 MB of RAM, 2x1 GB and 2x2 GB Hard Disk and 40 .Mflops 

( scalar ) peak performance. 

The performance of multigrid technique is assessed by cr*T::> serving the two parameters, 
namely. Central Possessing Unit ( CPU ) time and the Worlcx 3 . 33 -it. The Workunit is defined 
as the one relaxation sweep on the finest grid level. Table 5. - shows the CP U times and 
the Workunits for one gridsize ( 314 x 66 ) and different ^f^^ynolds numbers for steady 
flow. In the present computation, we have recorded the tota,l xxnmber of relaxation sweeps 
( total Workunits ) in the finest grid level for multigrid a.s -vv^ll as for the single-grid till 
the solution is converged. Comparison of both the treatme: 3 ^hs shows that the multigrid 
requires very less number of relaxation sweeps ( on finest gri<l level ) compared to single- 
grid. The CPU time for the multigrid varies in the range of 50 - 65 % of single-grid time. 
The saving in computational time due to the multigrid is ci,l>o'u.t 35 - 50 % in comparison 
to single-grid for different Reynolds numbers. 

The reduction of work obtained by the multigrid treat depends on the amount 
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RtB 

~ pUavR j P 


Workunits 
to get converged 
steady solution 

Amount of 
work in 
relation to 
single grid 

33 

Multigrid 

4239 

37.77 % 

Single- grid 

11223 

50 

Multigrid 

4974 

39.65 % 

Single-grid 

12545 

70 

Multigrid 

5712 

40.08 % 

Single- grid 

14253 

100 

Multigrid 

7317 

41.67 % 

Single-grid 

17560 


Table 5.4: Comparison of Workunits between the multigrid and single-grid for the three- 
dimensional problem 

of work for the corresponding single-grid computation of the continuity equation, since the 
amount of work for the momentum equation remains unchanged. The higher the fraction 
of the work for the continuity equation in comparison to the whole work, the higher is 
the reduction of the work achieved by the multigrid treatment. Our computational tests 
with multigrid show a reduction of about 80 - 90 % of work for a blockage ratio of 0.25 
depending on the Reynolds numbers. 

For the unsteady flows, we observed that for the multigrid calculation, the period- 
icity is set in the flow field at the expense of much less CPU time as compared to the 
single-grid calculations. 

5.2.2 Three-Dimensional Case 

In the multigrid treatment, the finest grid size is as the single-grid ( 82 x 26 x 26 ) 
for three-dimensional flows. Three grid levels have been deployed for computation. In 
Chapter-4, it has already been mentioned that for adaptive cycling, the switching factor 
6y is chosen as 0.6. This particular value can be varied from problem to problem. For 
the two-dimensional case, stated in the previous section, Sy equal to . 0.6 gives good 
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convergence rate but for the present configuration of interest, he, for a three-dimensional 
channel-flow, the optimum value is found to be 0.3. The other treatments concerning the 
operations of the multigrid technique have already stated in Chapter-4. 

The performance of the multigrid technique is not found very exciting in the three- 
dimensional case. The performance of the multigrid treatment basically depends on 
the type of errors that to be smoothened during relaxation. For the three-dimensional 
problem, it has been observed that the low frequency errors are smoothened out very 
quickly and only for the first two or three time steps it takes large number of relaxation 
sweeps on pressure-velocity equation to make the flow field divergence free. After the 
solution smoothened out in two to three time steps, successive time steps need only five 
to seven relaxation sweeps to make al the cells divergence free. In the multigrid method, 
every V-cycle needs atleast three sweeps in the finest level. However, some enhancement 
in computation speed has been observed for the multigrid calculations. Table 5.4 shows 
the comparison of Workunits between the multigrid and the single-grid calculations for 
the steady flow cases. 



Chapter 6 


Concluding Remarks 


Full Navier-Stokes equations have been solved for two and three-dimensional flows in 
a channel with a built-in square cylinder. Vortex shedding behind the square cylinder 
induces periodicity in the flow field. Details of the phenomenon are simulated through 
numerical flow visualization. Computations for three-dimensional flow show formation 
of horseshoe vortices in the wake of the square obstacle. The periodicity of the flow is 
characterized by the nondimensional frequency known as Strouhal number. 

A V-cycle multigrid technique has been suitably adapted to accelerate the speed 
of convergence of the solution. The multigrid technique improves the computation time 
significantly. However, for two-dimensional flows the multigrid technique appears to be 
more effective as compared to three-dimensional cases. For three-dimensional calcula- 
tions, the convergence becomes slow after the removal of high frequency errors. 
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